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Abstract 



Su-Doku, a popular combinatorial puzzle, provides an excellent test- 
bench for heuristic explorations. Several interesting questions arise from 
its deceptively simple set of rules. How many distinct Su-Doku grids are 
there? How to find a solution to a Su-Doku puzzle? Is there a unique 
solution to a given Su-Doku puzzle? What is a good estimation of a 
puzzle's difficulty? What is the minimum puzzle size (the number of 
"givens")? 

This paper explores how these questions are related to the well-known 
alldifferent constraint which emerges in a wide variety of Constraint Sat- 
isfaction Problems (CSP) and compares various algorithmic approaches 
based on different formulations of Su-Doku. 



Su-Doku grids and puzzles. Su-Doku is a well-known logic-based number 
placement puzzle. The objective is to fill a 9x9 square grid so that each line, 
each column or file, and each of the nine 3x3 blocks contains exclusively the 
digits 1 to 9, only once each. A puzzle is a partially completed grid. 

This definition is readily generalized to larger grids. Let us consider M n = 
{1 . . . n} the set of digits 1 to n, a Su-Doku of size n is a n 2 x n 2 grid which is 
to be filled so that each of the n 2 lines, each of the n 2 columns and each of the 
the n 2 blocks contains the digits 1 to n only one time each. 

The deceptively simple definition actually hides large combinatorial prob- 
lems, even for small values of n. This is, for instance, reflected in the number 
of such Su-Doku grids as calculated by Felgenhauer and Jarvis [6]: 



1 Su-Doku as a CSP 



Size 



Number 

1 

288 

6,670,903,752,021,072,936,960 



1 

2 
3 



Table 1. 



This is sequence A107739 in the Online Encyclopedia of Integer Sequences. 
The number of grids in the familiar size-3 Su-Doku is already mind-boggling! 



Constraint Satisfaction Problem. The rules of Su-Doku can be cast in 
terms of constraints that a solution should comply with to be valid. A constraint 
satisfaction problem (CSP) is a triple (X, D,C), where X is a sequence of n 
variables x\, X2, ■ ■ ■ x n , D is a sequence of n finite domains D\, D2, ■ ■ ■ D n , where 
Di is the set of possible values for variable Xi , and C a finite set of constraints 
between variables. 

A constraint C on a subsequence of variables Xi t , Xi 2 , . . . xi m is simply a 
subset of the cartesian product D il x D i2 . . . D im , which expresses the allowed 
combination of variable values. 

Constraints are often expressed as equations that variables must satisfy. 

General CSP Expression for Su-Doku. As will be seen in later sections, 
there are several ways to express a Su-Doku puzzle as a concrete CSP. In 
somewhat abstract terms, the CSP formulation would state four groups of con- 
straints: 

i. Each of the n 4 cells contains a unique value from the set M n i 

ii. Each of the n 2 lines of the grid contains all values from the set M n i 

iii. Each of the n 2 files of the grid contains all values from the set M n i 

iv. Each of the n 2 blocks of the grid contains all values from the set M n i 

where each of n 4 variables x\, . . . , x n 4 would represent a grid cell. 

The initial "givens" of a Su-Doku puzzle restrict the actual sets from which 
the above constraints allow values to be assigned to these variables. In the 
conventional Su-Doku puzzle, there are 81 variables and 81 + 3x9 = 108 such 
general constraints. 

1.1 Constraint propagation v. Search 

Alternating propagation and search. When nothing much is known about 
the constraints of a CSP, the general procedure to find a solution, or to show 
that none exists, alternates propagation and search. Each constraint is asso- 
ciated with a propagation procedure which tries to shrink the domains of the 
constraint variables by removing values that are certainly not part of a solution. 
This procedure uses first the locally available information, i.e. information pro- 
vided by the constraint itself. When such a value is excluded from a variable's 
domain, this information is also provided to all other propagation procedures 
associated with constraints that share the same variable. This externally pro- 
vided information might then trigger a propagation procedure of a constraint 
which was locally consistent, which in turns might exclude other values, provid- 
ing again this new information to other constraint procedures. The cascade of 
propagation procedure calls stops, however, after some time, in a state of global 
stability where all of the CSP constraints are locally consistent together. 

The globally stable state may or may not be a solution (or a failure) to the 
CSP. If not, a search procedure splits the current problem into at least two 
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CSPs, usually by assigning a value to a variable from its current domain. More 
generally this is done by dividing the current Di domain of a variable Xi into 
k disjoint subsets (k ^ 2) the union of which is Di and considering the k CSP 
subproblems. The propagation phase is repeated for each subproblem, thus 
recursively creating a search tree with the original problem at its root and the 
solutions or failures at its leaves. It is indeed a search procedure since the choice 
of the domain to split, and hence of which variable to propagate, and how to 
split it may be based on heuristics or rules of thumb specific to the problem at 
hand. 

Human Problem-Solving in Su-Doku puzzles. It is generally acknowl- 
edged that the alternating propagation and search procedures are effectively 
used by human solvers when tackling Su-Doku puzzles. All tutorials, books, 
Web sites and other broadly available Su-Doku instructional material start by 
highlighting that the first steps to be taken towards a solution involve propaga- 
tion of the givens, simply by ticking off for each of them its value from the other 
cells in the same line, file and block. Should this leave exactly one possible 
value for another cell, mark it and propagate in turn to this cell's line, file and 
block until no more propagation is possible. 

At this step, there are litterally tens of heuristics with inventive names such 
as "Fish" , "XY-Wing" or "XY-Chain" , ranging from the simple check to the 
complex pattern, to choose from in order to perform the search phase. Numer- 
ous Su-Doku Web sites offer catalogs of patterns to look for in a partially solved 
puzzles, and each human solver develops his or her own individual catalog as 
experience grows. 

There is no recognized standard, however, for evaluating the complexity of 
such heuristics. Hence there is no simple way of comparing Su-Doku puzzles 
difficulties, as some search heuristics may be considered simple by some and 
complex by others. While the number of givens, which are in fact starting 
points for the propagation procedures, is certainly an indication of the difficulty 
level, it is not a complete indicator of the actual complexity. Even so-called 
minimal puzzles where for n — 3 there are 17 givens, might prove easy to solve 
if the initial propagation goes far towards the solution, leaving only few empty 
cells for search. (The fact that a n — 3 Su-Doku puzzle needs at least 17 
givens in order to have a unique solution is still, at the time of this writing, a 
conjecture.) 

1.2 A Simple Propagation-Search Algorithm 

In this section we present a very simple implementation of the alternated propagation- 
search phases to solve Su-Doku puzzles as CSP. This is for illustrative purpose 
and by no means the only way to implement propagation and search, or to strike 
a balance between propagation and search in CSP solutions. Some of the ideas 
here are inspired by [15] . and, for lack of a better name, we simply call this 
algorithm the PS-1-2 algorithm. 
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Propagation. With each cell in the grid, the algorithm maintains an array 
of the valid values which can be used for this cell, its so-called domain that the 
propagation phase seeks to reduce as much as possible provided the constraints. 

Initially for a n order Su-Doku puzzle, all domains Dij are the same set 
M n 2 of the first n 2 integers. 

Propagation resolves into iterating four separate steps: 

— A value v is assigned to a cell, thus reducing Dij to {v} 

— The newly assigned value v is deleted from the domains of cells lying in the 
same line, file and block than the initial cell: -Di,fc, D^j and fs(i,j)=s(p,g) 
are therefore reduced to Di t k\{v}, Dkj\{v} and D B ^ij^ =B ^ pq ^\{v} re ~ 
spcctively. (Here i,j range over 1 . . .n 2 , and so do k,p and q; the some- 
what informal notation B(i,j) denotes the block, in the range l...n 2 , 
containing the i,j cell.) 

— Another reduction step is taken, leveraging the duality of constraints in 
the particular Su-Doku CSP formulation. As each line, file or block may 
only contain one occurence of any number 1 . . . n 2 , each of the number has 
to be assigned to a unique cell within a line, file or block. So if a value 
v appears only once in all Dij of a given line, file or block, the process 
reduces this singled out Dij to {v}. 

— After this reduction step, some domains Dij may happen to be empty, in 
which case the puzzle has no solution - we'll say that the propagation is 
blocked - or to be reduced to a singleton {w}. If all domains are singletons, 
the puzzle is solved. If only some of the domains are singletons, the 
algorithm reiterates these steps assigning the singletons' values to the 
corresponding cells. 

The iteration is stopped when no further reduction happens in step 4 of the 
above propagation process. Reductions are done in any order as it does not 
impact the final result after the system reaches a quiescent state. The "1" 
in the algorithm name comes from the choice of reducing domains on a single 
constraint type (and its dual): the unicity of values for CSP variables. 

Data representation. In order to lower the computation costs, the domains 
for each of the n 2 variables representing the puzzle cells are implemented as 
packed arrays in C. Reduction then becomes a logical operation on a bit ar- 
rays. Step 3 of the previous propagation process requires the domains to be 
transposed: for each line, file and block, n 2 new bit arrays are computed, the 
i-th of which is made of bits i of the n 2 domain bit arrays. 

Binary Search. After the propagation phase we may end in a blocked, solved 
or yet indeterminate state. 

Again we choose here a simple search procedure for the next steps towards 
a solution in the indeterminate case. Namely we look for domains reduced to 
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simple pairs {i>, w} and we operate a binary depth-first search, first reducing 
to {v} and then, if no solution is found after recursive propagation and search, 
reducing to {w} and propagating/searching again. 

Note that if at the end of each of the propagation phase we do indeed find 
pair domains, this process will certainly terminate either finding a solution or 
proving, after enumeration, that no solution exist for this puzzle. 

The "2" in the name of the algorithm derives from the fact that we only 
consider pair domains in the search phase. 

The PS-T2 algorithm is only one in a scope of algorithms which we will in- 
vestigate further in the following sections and which were developed to solve 
a particular type of constraint, called alldifferent in the CSP literature. It so 
happens that the CSP formulation of Su-Doku uses only inter-related alldiffer- 
ent constraints thus offering a perfect case study for combinatorial analyses of 
the various approaches to the general solution of alldifferent constraints. Even 
more so, we will see that Su-Doku constraints are a special form of alldifferent 
constraints, i.e permutation constraints, which specific properties will suggest 
a completely different approach to Su-Doku puzzle representation and solution 
explored in the second section of this paper. 

1.3 Experiments and Results 

A C Implementation. The PS-T2 algorithm was implemented in C under 
Cygwin for experimentation purposes. The core of the implementation is ar- 
ticulated around two functions: a propagation and reduction function called 
solveStep, and a recursive depth-first search function called pairReduce. 

int solveStep ( int main_step ){ 
int step, i, flag, main_flag = 1; 
while ( main_flag ){ 

// Local rules and propagation main loop 

flag =1; 

step =1; 

// 1. Propagate givens 
while ( flag ){ 

flag = propagateQ; 

step++; 

} 

main_step += step; 

// 2. Reduces lines, cols and blocks 

// Rem.: flag is at this point, L is n*n 

for( i = 0; i < L; i++ ) flag += reduceLine( i ); 

for( i = 0; i < L; i++ ) flag += reduceColumn( i ); 

for( i = 0; i < L; i++ ) flag += reduceBlock( i ); 

main_step += flag; 
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main_flag = (flag > ) ? 1 : 0; 
} 

return main_step; 

} 

The previous code fragment details the solveStep function which propagates 
assignments of values to cells by calling the (not-represented) propagate func- 
tion, which in turn operates on the domain bit array representations, deleting 
the assigned values from other cells' domains in each relevant line, file and 
block. This is in fact step 2 of the PS- 1-2 algorithm as described in the pre- 
vious section. Then the dual step in domain reduction is taken by calling 
the (not-represented) reduceLine, reduceColumn and reduceBlock functions 
which handle the transposition and reduction in step 3 of the PS-1-2 algorithm. 

This function exits when no domain can be further reduced to a singleton 
through the iteration of the basic propagate and reduce operations. 

In addition the function maintains various counters, namely step and main_step, 
for simple statistics. 

The depth-first search function is straightforward: 

int pairReduce( int step ){ 
packPtr p; 

struct pack keep[S]; 

if ( 1 == solvedpO ) return step; 

if( 1 == blockedpO ) return step; 

S_ReduceSteps += 1; 

// Find a pair domain p 

p = nextPair () ; 

if ( (packPtr)O == p ){ 

printf( "No pair left. S=°/„d, B=7.d 
", solvedpO, blockedpO ); } 
else{ 

int hi, lo; 

// Store current state of search as an array of domain bit arrays 
packcpy( keep, cell ); 

// Extract low and high values in pair p 

lo = getPack( p ) ; 

// Delete low value from domain p 

subPack( p, lo ) ; 

hi = getPack( p ) ; 

// And propagate to other domains 

step = solveStep ( step ); 

if ( 1 == solvedpO ) return step; 

step = pairReduce( step ); 

if ( 1 == solvedpO ) return step; 
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//No solution reached, restore search state 

packcpy( cell, keep ); 

// Now delete high value from domain p 

subPack( p, hi ) ; 

// And propagate to other domains 

step = solveStep ( step ); 

if ( 1 == solvedpO ) return step; 

if ( == blockedpO ){ 

step = pairReduce( step ); 

} 

} 

return step; 

> 

This code fragment illustrates the search procedure. If pairReduce is entered in 
a solved or blocked state it returns immediately. If entered in an indeterminate 
state, it first finds a pair domain, by calling the (non-represented) nextPair 
function. (In the implementation this function does a simple but costly linear 
search on the array of all domain bit arrays.) If it fails to find such a pair 
domain it simply stops, although in all of the test puzzles this never happened. 

When it succeeds, however, the function backs up the current search state, 
here an array of domain bit arrays representing the remaining possible values 
for each cell in the puzzle, assigns first the highest value of the pair domain 
to the cell and propagates this assignment by calling the previously mentioned 
solveStep. At this point, the puzzle is either solved and pairReduce returns 
(recursively up to the first caller in fact) , or in a blocked or indeterminate state. 
In the latter cases, and in standard depth-first fashion, we search another pair 
domain by calling recursively pairReduce. 

Note that if the state is blocked, this recursive call returns immediately. 
When it does not and we still have no solution on the first branch of our binary 
search, we try the other one. The function duly restores the backed up state of 
search and this time assigns the lowest value of the pair domain to the cell and 
propagates to other domains, again calling solveStep. 

Various counters are also updated for statistical purposes. In order to 
obtain a running solver program, these functions are wrapped into a main func- 
tion which initializes all domains to the same all-one bit arrays, reads the puz- 
zle in from a file, executes the initial propagation of the "givens" and calls 
pairReduce (0) . 

An Example Run. We ran PS-1-2 on some of the minimal puzzle instances 
as collected by Gordon Royle [TH] who maintains a catalog of order 3 Su-Dokus 
puzzles with only 17 "givens" . 

On another example, the following puzzle, which is not part of this "minimal 
puzzles" set: 
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Example 1 Puzzle 
.125.487. 



75 23 

. .41.87. . 
.2. .5. .4. 
. .34.95. . 
48 17 



.357.169. 

is solved with only 77 propagation and 11 search operations by PS- 1-2, as 
detailed in the following tabular trace: 



Example 2 
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12 
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23 
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2 









26 


2 






h 


26 
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3 
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3 









29 


3 
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29 


4 


5 






34 


4 




1 




35 
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3 






38 


4 









38 


4 
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38 
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3 






41 
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41 


5 
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41 
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5 






46 


6 




1 




47 


6 
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50 
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50 


6 
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50 


6 



8 



3 






53 
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53 
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h 


53 
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3 






56 


7 




1 
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60 


7 









60 


7 
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60 


8 


5 






65 


8 









65 
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65 
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68 


9 









68 


9 






h 


68 


10 


4 
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10 
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10 
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11 


5 






77 


11 









77 


11 



Grid: solved 1, blocked 0; in 77 operations 
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1 


7 I 
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7 
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2 
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8 I 


1 2 
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5 


8 


1 


6 


3 


9 


4 I 



The process called the search procedure 11 times, when the propagation/reduction 
operations reach quiescence as indicated by a in the Red(uctions) column. The 
Srch column indicates whether the h(igh) or l(ow) value of the pair searched 
is used for the next propagation phase. In the particular instance, backtrack 
occurred only once at the sixth pair search: both high and low value were 
propagated to find the solution. 

Conclusions. The canonical procedure to solve CSP-formulated problems al- 
ternates a propagation phase, where data is used to reduce domains of the 
variables as far as possible, also known as filtering, with a search phase, a back- 
track procedure which explores incremental steps towards a solution. There 
is ample room for variability in this framework both in the balance between 
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propagation and search, and within each phase in the criteria used in filtering 
and in search. 

In the case of Su-Doku puzzles, we have presented a naive algorithm, PS-1-2, 
which only filters on unicity of the variable value and of this value per group 
(line, file or block) in the propagation phase, and only uses binary search in the 
alternating search phase. Although there should be pathological cases where the 
binary search phase might fail, the PS-1-2 algorithm was successful at solving 
quickly all the puzzles we submitted, including so-called minimal puzzles. 

1.4 CSP/SAT/LP formulations: the alldifferent constraint 

While several ad hoc CSP solving procedures may be designed for Su-Doku 
puzzles, its constraints generally fall under a now well documented class of 
constraints for which efficient filtering procedures have been published in the 
literature and are embedded in several tools, commercial and otherwise. The 
pattern appearing in all the constraints in the above CSP formulation of Su- 
Doku directly relates to one of the latter, the alldifferent constraint [2T] , 

A CSP alldifferent expression. We will rephrase the CSP expression in 
terms of this well studied alldifferent constraint. Let us consider the n-sized 
Su-Doku puzzle and introduce the n 4 variables xi t i . . .x n 2 n 2 representing the 
cells in a grid where, by convention, Xij, is the value to be assigned to the cell 
in line i and file j. All variables have the same domain, taking their values in 
M„2. The CSP expression of the problem is to find a unique value for each 
variable satisfying the following set of alldifferent constraints: 

• Vi S {1, . . . n 2 } alldifferent(a;i ! i, . . . x n 2 

• Vj £ {1, . . . n 2 } alldifferent (xj.i, . . . x^ n i) 

• V& E {1, . . • n 2 } alldifferent (xbj , . . . Xf, 2 ) where the bk are the pair of indices 
of variables representing cells in the same block 

SAT formulations. A given alldifferent constraint naturally translates into 
a set of simpler binary constraints on its variables, the naive translation. In 
such naive translation the alldifferent constraint is expressed as a conjunction 
of disjunctive clauses involving at most two variables and the values of the 
variable domains. 

For instance in the size 2 Su-Doku puzzles, an alldifferent constraint on four 
variables as in line, file and block constraints becomes: 
alldifferent (x\, x 2 , X3, £4) ^ 

(xx = 1 V xi = 2 V xi = 3 V x 1 = 4) A (x 2 = 1 V x 2 = 2 V x 2 = 3 V x 2 = 4) A 
(x 3 = 1 V x 3 = 2 V 23 = 3 V x 3 = 4) A (x 4 = 1 V x A = 2 V x A = 3 V x 4 = 4) 
a CNF formula with 10 disjunctive clauses, 6 of which are binary and 4 
of which unary. Generally speaking the size n alldifferent constraint naively 
translates to a CNF formula with 77,(77, + l)/2 disjunctive clauses. 
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The naive translation is actually enough to be fed to standard SAT solvers 
such as maxsatz [10] and minisat [5], for instance. 

LP formulations. Logic programming tools can also directly use the above 
SAT expressions. In this section we investigate the use of the CLP(FD), or Con- 
straint Logic Programming for Finite Domains, extension to the GNU-Prolog 
implementation [1] in Su-Doku puzzle experiments. 

Expanding on the above analysis, a naive implementation of a single alld- 
ifferent constraint simply translates it into a corresponding set of binary con- 
straints, each one of which stating that a given variable is different from the 
other. Focusing on 4x4 Su-Doku grids, i.e. n = 2, for instance, such a naive 
implementation of the single alldifferent constraint on four variables would then 
be as follows: 

naive_all_different(X,Y,Z,ZO) :- 

X \= Y, X \= Z, X \= Z0, Y \= Z, Y \= ZO, Z \= ZO . 

which states that each variable should have a distinct value from the other 
three variables in the group. In order to complete the size 2 Su-Doku grid 
enumeration, a definition predicate, assign, is created to define the (unique) 
domain of all variables: 

assign(l) . 
assign(2) . 
assign(3) . 
assign(4) . 

It is a predicate which is true for each of the four admissible values of the 
cells in a size 2 Su-Doku puzzle. To complete the GNU-Prolog program, we 
use these two predicates to express all the constraints of the 4x4 grid: 

naive_puzzle( A00, A01, A10, All, BOO, B01, BIO, Bll, 
COO, C01, CIO, Cll, D00, D01, D10, Dll ) :- 
system_time (TO) , 
cpu_time (T10) , 
real_time(T20) , 



assign( 


A00 


), 


assign( 


A01 


) 


assign( 


A10 


) 


assign( 


All 


) 


assign( 


BOO 


) 


assign( 


B01 


) 


assign( 


BIO 


) 


assign( 


Bll 


) 


assign( 


COO 


) 


assign( 


C01 


) 


assign( 


CIO 


) 
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assign( Cll ) 

assign( DOO ) 

assign( D01 ) 

assign( D10 ) 

assign( Dll ) 



naive_all 
naive_all 
naive_all 
naive_all 

naive_all 
naive_all 
naive_all 
naive_all 



different 
different 
different 
different 

different 
different 
different 
different 



(A00, A01, 
(BOO, BOl, 
(COO, C01, 
(DOO, D01, 

(A00, A01, 
(A10, All, 
(COO, C01, 
(CIO, Cll, 



A10, All) 

BIO, BID 

CIO, Cll) 

D10, Dll) 

BOO, BOl ) 

BIO, Bll ) 

DOO, DOl ) 

DIO, Dll ) 



naive. 


_all_dif f erent (A00 , 


A10, 


COO, 


CIO 


) 


naive. 


_all_dif f erent (A01 , 


All, 


C01, 


Cll 


) 


naive. 


_all_dif f erent (BOO , 


BIO, 


DOO, 


DIO 


) 


naive. 


_all_dif f erent (BOl , 


Bll, 


DOl, 


Dll 


) 



system_time (T) , 
cpu_time (Tl) , 
real_time (T2) 

write( 'time TO: '), write(TO), write (', time T: ' ),write(T), nl , 
write ( 'time TO: '), write(TlO), write (', time Tl: ' ),write(Tl), nl, 
write ( 'time TO: '), write(T20), write (', time T2 : ' ),write(T2), nl . 

The predicates assigns unique values to the 16 variables representing the 
corresponding cells of the Su-Doku grid or puzzle, and expresses the alldifferent 
constraint on each of the 4 lines, files and blocks. It also keeps track of various 
execution times for instrumentation purposes. 

Execution times are unsurprisingly very long, even for the size 2 Su-Doku 
grids as the naive implementation only uses extremely local constraints within 
one line, file or block and ignores global constraints. 

Fortunately, GNU-Prolog bundles a constraint logic programming extension 
for finite domains which incorporates the latest filtering algorithms for a wide 
variety of constraints used in industry problems. The CLP(FD) extension uses 
specific filtering techniques for the alldifferent constraint, the theoretical basis 
of which will be explored in the next sections, leading to a much more efficient 
implementation of enumeration and solving tasks. 

More specifically, the CLP(FD) extension offers a simple set of built-in pred- 
icates such as: fd_domain to define variables' domains, fd_all_diff erent to 
state a single alldifferent constraint (of any arity), and fd_labeling to trigger 
search according to a choice of filtering methods. The CLP(FD) is fully de- 
scribed in Diaz's Thesis [3]. The size 2 Su-Doku implementation now becomes: 

puzzle( A00, A01, A10, All, BOO, BOl, BIO, Bll, 
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COO, COl, 


CIO, Cll, 


D 


f A 

id. 


.domain 


ADD 


i , 


4 


■PA 


.domain 


' AA1 
AU 1 


i , 


4 


-f A 


.domain 


' a 1 n 

A1U 


i , 


4 


■PA 


.domain 


' A 1 1 

Al 1 


i , 


4 


f d. 


.domain 


; boo 


i , 


4 


fd. 


.domain 


: boi 


i, 


4 


fd. 


.domain 


; bio 


i, 


4 


fd. 


.domain 


; bii 


i, 


4 


fd. 


.domain 


; coo 


i, 


4 


fd. 


.domain 


; coi 


i, 


4 


fd. 


.domain 


; cio 


i, 


4 


fd. 


.domain 


; en 


i, 


4 


fd. 


.domain 


; doo 


i, 


4 


fd. 


.domain 


; doi 


i, 


4 


fd. 


.domain 


: dio 


i, 


4 


fd. 


.domain 


: dii 


i, 


4 



DOI, DIO, DID 



fd_all_different([AOO, A01, A10, All]) 

f d_all_dif f erent ( [BOO , BOI, BIO, Bll]) 

fd_all_different([COO, COl, CIO, Cll]) 

fd_all_different([DOO, DOI, DIO, Dll]) 

fd_all_different([AOO, A01, BOO, BOI ]) 

fd_all_different([A10, All, BIO, Bll ]) 

fd_all_different([COO, COl, DOO, DOI ]) 

fd_all_different([C10, Cll, DIO, Dll ]) 



diff erent ( [A00, 


A10, 


COO, 


CIO 


]) 


different ( [A01, 


All, 


COl, 


Cll 


]) 


diff erent ( [BOO, 


BIO, 


DOO, 


DIO 


]) 


different ( [BOI, 


Bll, 


DOI, 


Dll 


]) 


time (TO) , 











cpu_time(T10) , 
real_time(T20) 

fd_labeling([AOO, A01, A10, All 
BOO, BOI, BIO, Bll 
COO, COl, CIO, Cll 

DOO, DOI, DIO, Dll] , [variable_method(most_constrained)] ) , 
system_time (T) , 
cpu_time (Tl) , 
real_time(T2) 

write( 'time TO: '), write(TO), write(', time T: ' ),write(T), nl, 
write( 'time TO: '), write(TlO), write (', time Tl: ' ),write(Tl), nl 
write( 'time TO: '), write(T20), writeC, time T2 : ' ),write(T2), nl . 
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Basically this new puzzle predicate defines the 16 domains for the 16 variables 
with values in the same M4 set, states the three group of alldifferent constraints 
and finally searches for a proper labeling. As before, several ancillary predicates 
have been added for instrumentation purposes. 

On the simple enumeration task, the first solution is almost instantly com- 
puted on a standard Intel machine running Windows XP: 

puzzle( A00, A01, A10, All, BOO, B01, BIO, Bll, 
COO, C01, CIO, Cll, D00, D01, D10, Dll) . 
time TO: 296, time T: 296 
time TO: 1609, time Tl: 1609 
time TO: 155875, time T2 : 155875 



A00 = 1 
A01 = 2 
A10 = 3 
All = 4 
BOO = 3 
B01 = 4 
BIO = 1 
Bll = 2 
COO = 2 
C01 = 1 
CIO = 4 
Cll = 3 
D00 = 4 
D01 = 3 
D10 = 2 
Dll = 1 ? ; 

time TO: 296, time T: 312 
time TO: 1609, time Tl: 1625 
time TO: 155875, time T2 : 158472 



A00 = 1 

A01 = 2 

A10 = 3 

All = 4 

BOO = 3 

B01 = 4 

BIO = 1 

Bll = 2 

COO = 2 

C01 = 3 

CIO = 4 

Cll = 1 

D00 = 4 
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D01 = 1 
DIO = 2 
Dll = 3 ? ; 

time TO: 296, time T: 328 
time TO: 1609, time Tl : 1641 
time TO: 155875, time T2 : 222624 

A00 = 1 
A01 = 2 
A10 = 3 
All = 4 
BOO = 3 
B01 = 4 
BIO = 1 
Bll = 2 
COO = 4 
C01 = 1 
CIO = 2 
Cll = 3 
D00 = 2 
D01 = 3 
DIO = 4 
Dll = 1 ? ; 

time TO: 296, time T: 343 
time TO: 1609, time Tl : 1656 
time TO: 155875, time T2 : 228535 

and the other 288 solutions are quickly printed out for a total execution time 
of 2,797 milliseconds. 

1.5 Arc-Consistency and value graph 

In order to understand how to make the best use of the local information pro- 
vided by the constraints themselves, we introduce some definitions of local con- 
sistency. 

Definition 1 A constraint of arity m on the variables hyperarc 
consistent if all values of the variables are used in some solution to the con- 
straint, i.e. \/x ik \/v e D ik , 3(v h , . . . , v ikl , v ik+1 , . . . , v lm ) G D h x . . . D ik l x 
Di k+1 x . . . Di m such that (v^, . . . ,Vi k _ 1 ,v,Vi k+1 , . . . ,Vi m ) is a solution to the 
constraint. 

Definition 2 A constraint C is arc- consistent when C is of arity 2 (binary) 
and hyperarc consistent. 

A constraint has a solution if it can be made hyperarc consistent. Hyperarc 
consistency is the best possible pruning based on the local information provided 
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by the constraint. The naive implementation above translate the alldifferent 
constraint into a collection of binary constraints which are made arc-consistent 
by filtering the domains following the simple procedure which as soon as a 
domain is reduced to one value, removes this value from the domains of all 
other variables. 

The more efficient variants used, for instance, in the CLP(FD) package rely 
on a completely different approach based on results from graph theory. The 
correspondence with graph theory was used by Regin [T5] to create a filtering 
algorithm from matching theory. We introduce the notion of bipartite graph. 

Definition 3 Bipartite Graph. A graph G consists of a finite, non-empty set 
of elements V called nodes, or vertices, and a set of unordered pair of nodes E 
called edges. If V can be partitioned into two disjoint, non-empty sets X and 
Y such that all edges in E join a node in X to a node in Y, G is called bipartite 
with partition (X,Y); we also write G = (X,Y,E). 

The definition directly applies to the alldifferent constraint say, on variable 
set X = {x\, . . . x n } and domains D\ . . . D n , in that it specifies its value graph. 

Definition 4 Value Graph. Given an alldifferent constraint, the bipartite graph 
G = (X, (J Di , E) where (xi, d) G E iff d G Di is called the value graph of the 
constraint. 

The value graph has an edge from each variable in the constraint to each of 
its domain value. Solving such a constraint becomes a problem of maximum 
matching in the corresponding value graph. 

Definition 5 Maximum Matching. A subset of edges in a graph G is a match- 
ing if no two edges have a vertex in common. A matching of maximum cardi- 
nality is called a maximum matching. A matching covers a set of vertices X isf 
every node in X is an endpoint of an edge in the matching. 

The link between matching theory and hyperarc consistency established by 
Regin is as follows. 

Proposition 1 The constraint alldifferent on variable set X is hyperarc consis- 
tent iff every edge in its value graph belongs to a matching that covers X in the 
value graph. 

Hyperarc and arc-consistency algorithms are around 0(dn 15 ) where d is the 
maximum cardinality of the domains and n the number of variables. 

Matching Theory. Obviously if there is a complete matching from X to Y 
in (X,Y,E) then for every S C X there are at least \S\ vertices of Y adjacent 
to a vertex in S. That this necessary condition is also sufficient is usually 
called Halls' theorem. This fundamental result was proved by Hall in 1935, 
but an equivalent form of it had been proved by Konig and Egervary in 1931; 
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both versions, however, follow from Menger's theorem from 1927. We refer to 
Bollobas for demonstrations and historical remarks pQ. 

These results will also be the origin of yet another approach to solving Su- 
Doku puzzles and enumerating grids, as a complete matching is also called a 
set of distinct representatives, from which are derived new expressions of the 
Su-Doku problems in terms of exact covers or dually exact hitting set problems. 
This new formulation will suggest a different algorithmic approach which is 
explored in Section 2. 

Regin's algorithm relies on the fact that if we know only one arbitrary maxi- 
mum matching, we can efficiently compute if an edge of the value graph belongs 
to some matching of the same maximum size without having to explore all such 
matchings. 

1.6 Bounds- and range-consistency filtering 

Hall's theorem from matching theory may also be used in relation to weaker 
forms of consistency called bounds-consistency and range-consistency. 

Definition 6 Bounds Consistency. A constraint of arity m where no domain 
Di is empty, is called bounds consistent iff for each variable and each value in 
the range bounded by min(Z?i) and max(Di), there exist values in the respective 
ranges bounded by the other domains minimum and maximum values such that 
together with the latter Di value they constitute a solution to the constraint. 

Definition 7 Range Consistency. A constraint of arity m where no domain 
Di is empty, is called bounds consistent iff for each variable and each value in 
Di, there exist values in the respective ranges bounded by the other domains 
minimum and maximum values such that together with the latter Di value they 
constitute a solution to the constraint. 

Note that here all domains are supposed to be integer domains which can 
be (totally) ordered. In contrast to hyperarc and arc-consistency, bounds and 
range consistency look for values in the intervals defined by domains, rather 
than the domains themselves. As these intervals may be larger than the actual 
domains, these notions represent weaker form of consistencies than hyperarc 
and arc-consistency. Both may be considered as a relaxation of the hyperarc 
consistency. (In addition, bounds consistency may be regarded as a relaxation 
of range consistency itself.) 

Hall's theorem has been applied by Puget to create a bounds consistency 
algorithm for the alldifferent constraint [T7]- Given an interval / let us denote 
Kj the set of variables Xj such that Dj C /, i.e. the subset of variables which 
domains are included into the considered interval. We say that / is a Hall 
interval iff |/| = \Kj\. Puget 's result is as follows. 

Proposition 2 The constraint alldifferent on variables x%, . . . , x n where no do- 
main Di is empty is bounds consistent iff 
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i. for each interval I \Kj\ < \I\, 

it. for each Hall interval I , {min(Di), max(_Di)} n I = for all Xi $ Kj. 

This can be used to create an algorithm for bounds consistency on alldifferent 
constraints. We check every interval I with bounds ranging from the minimum 
of all domains to the maximum of all domains. When |/| < \Kj\ the constraint 
is inconsistent. And for each Hall interval, we remove the minimum and/or 
maximum until the intersection with I is empty. 

Example 3 Consider the following constraint: 

1. alldifferent (xi, X2, X3) 

2. D x = {1, 2}, D 2 = {1, 2}, and^g = {2, 3} 

The intervals we need to check are [1, 2], [1, 3], and[2, 3]. When I = [1, 2], the 
domains included in the interval are D\ and D2, hence Kj = {x\, X2} and since 
|/| = \Kj\ = 2, I is a Hall interval. We only have one variable not in Kj, 
namely £3, for which {min(Z?3), max(Z?3)} n I = {2}. The algorithm removes 
then 2 from D3 and the resulting system of sets {1, 2}, {1, 2}, {3} is now bounds 
consistent. The two solutions (1, 2, 3) and (2,1,3) are now a simple consequence 
of the reduction of D3. 

Faster implementation of bounds consistency have been designed since Puget's 
publication [TTJ [T3] . Leconte introduced an algorithm that achieves range con- 
sistency [5], also based on the Hall's theorem. In dual definitions from above, 
given a set of variables K let Ijc be the interval [min(£>/<-), max(Dx)] where 
Dk is the union of all variable domains; if is a Hall set iff \K\ = \Ik\- 

Proposition 3 The constraint alldifferent alldifferent on variables x±,...,x n 
where no domain Di is empty is range consistent iff for each Hall set K, Di fl 
I K — for all Xi K. 

An algorithm for range consistency can be derived from the previous propo- 
sition in a similar way to the derivation of the bounds consistency algorithm. 

1.7 Su-Doku as a CSP: a specific problem 

Enumerating Su-Doku grids or solving puzzles within reasonable space and time 
limits requires efficient algorithms for a single type of constraint, the alldifferent 
constraint. While general "propagate + search" constraint-solving algorithms 
will work well on Su-Doku puzzles, specifically tailored algorithm for the alldif- 
ferent constraint work still better. There are several versions of such algorithms, 
which embody different degrees of consistency and hence of performance in the 
puzzle task. Modern CSP, SAT and LP-solvers usually provide specific filtering 
procedures for the alldifferent constraint, which make them tools of choice for 
studying the complexity of the Su-Doku universe. 
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2 Su-Doku as an Exact Cover problem 



In fact, Su-Doku grids and puzzles involve an even more specific constraint 
than the alldiffcrent constraint. As is obvious from the previously mentioned 
SAT and CSP formulations, the alldiffcrent constraints in the Su-Doku problems 
involve variables having the same domain, namely M n i . In addition the number 
of variables in a given constraint is equal to the size of their shared domain. We 
call such special alldifferent constraint, permutations. 

In this section we explore a completely different approach to enumerating 
and solving grids and puzzles based on the precedent observation. Although 
the theoretical basis, going back to matching theory and Hall's results in graph 
theory, is the same, the resulting algorithms will significantly differ from the ones 
derived from consistency-checking filtering procedures explored in the previous 
section of this paper. 

2.1 The Exact Cover and Exact Hitting Set problems 

Definition 8 (Exact Cover) Given a family 21 = {A\, . . . ., A m \ of subsets of a 
set X, 21 is called an exact cover of X when Vx 6 X, 3i, 1 ^ i < m such that 
x E Ai 

This definition has also a dual formulation which describes an exact hitting 

set. 

Definition 9 (Distinct Representatives) Given a family 21 = {Ai, . . . ., A m } of 
subsets of a set X, a set of m distinct elements of X, one from each Ai is called 
a set of distinct representatives of 21, or a hitting set. 

If the set of distinct representatives is X , it is called an exact hitting set and 
21 is an exact cover of X . 

Finding exact hitting sets and enumerating exact hitting sets may be solved 
by backtracking algorithms. Their efficiency in this case relies on the fact that 
a representative has a unique value among the possible ones in the Ai subset 
and the filtering procedure, in the previous section sense of CSP solving, then 
reduces to the simple elimination of this value from all other domains. 

A permutation problem is a constraint satisfaction problem in which each 
decision variable takes an unique value, and there is the same number of values 
as variables. Hence any solution assigns a permutation of the values to the 
variables. There are to! such permutations for a constraint involving to vari- 
ables. The important feature of permutation CSPs is that we can transpose 
the roles of values and variables in representing the underlying problem to give 
a new dual model which is also a permutation problem. Each variable in the 
original (primal) problem becomes a value in the dual problem, and vice versa. 

An injection problem is a CSP in which each decision variable takes a unique 
value, and there are more values than variables. (Obviously if there are fewer 
values than variables, the problem is trivially unsatisfiable.) 
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The primal and the dual permutation problems are of course equivalent, 
but efficient algorithm can leverage this transposition by switching from one 
model to the other when appropriate [TJ. Actually the simple PS-1-2 algorithm 
described in the previous section indeed used the fact that the specific alldifferent 
constraints in Su-Doku problems are all permutation problems: the so-called 
reduce functions transpose values to variables, seeking to further reduce primal 
and dual domains when possible. 

In Su-Doku problems, each constraint is both a permutation and an exact 
hitting set problem as all variables have the same domain. 

Matrix representation of exact hitting set and exact cover problems. 

Both exact hitting set and exact cover problem can be represented as follows. 
Given a boolean matrix M with n rows and m columns, the problem is to find 
a subset 21 of the rows of M, such that each column j in M has exactly one row 
i e 21 such that My = 1. 

For example, consider the following matrix M with n — 6 and m = 4: 
/ 1 1 \ 
11 



1 
\ 1 1 / 

Possible exact hitting sets are {1,3,5}, {2,6} and {3,4}. 
In the case of permutation problems, the matrices involved are square n x n 
matrices. 

2.2 Enumeration 

The count of exact hitting sets is the number of solutions to the constraints used 
in Su-Doku formulations. Generally speaking, the number of exact hitting sets 
for permutation constraints, i.e. in which the number of values is the same as 
variables, is given by the permanent of the representation matrix [12j . 

Definition 10 (Permanent) If A is an n-square matrix then the permanent of 
A is defined by 

n 

per A = ^2 

ff£S„ z=l 

where the summation extends over S n , the symmetric group of degree n. 

The permanent is an appropriate invariant for matrices that arise in combi- 
natorial investigations where the problem is essentially unaltered by relabeling 
of the items under consideration, which is obviously the case in permutation 
problems. For example, the total number of derangements ("le probleme des 
rencontres") of n distinct items is given by per(J — I n ) where J is the n-square 
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matrix with every entry equal to 1, and is the n-square identity matrix [16] . 
(This is series A000166 in the Online Encyclopedia of Integer Sequences.) 

Even though the permanent looks superficially likes the more familiar de- 
terminant (without the alternating ± signs), Polya observed that no uniform 
affixing of ± signs to the elements of the matrix can convert the permanent 
into the determinant, for n > 2. The apparent simplification of definition from 
the determinant results counter-intuitively in tremendous complications in the 
evaluation of permanents. In particular, and in contrast to the determinant, 
the permanent is not well-behaved under permutation of rows and columns of 
the matrix; it is, however, multilinear like the determinant. 

Van der Waerden's Conjecture. Bounds for the permanent have been 
found, however difficult its exact computation turns out to be. Given a n-square 
matrix A, the z-th row sum of A is defined by 

71 

T{ = (H,j 
3 = 1 

and similarly the i-th column sum of A is defined by 

n 

C i = a i,3 
3 = 1 

With these, we introduce a doubly stochastic matrix with this definition. 

Definition 11 (Doubly Stochastic Matrix) Let A = (ai.j) be a n-square matrix, 
then A is doubly stochastic if 

^ a-i.j ^ 1 for 1 ^ i,j ^ n 

Ti = 1 for 1 ^ i ^ n 

Ci = 1 for 1 ^ i ^ ri 

Note that the representation matrix of an exact hitting set (or exact cover 
problem) is amenable to a doubly stochastic matrix, in the case of permutation, 
by replacing each entry equal to 1 with 1/n. 

Van der Waerden made a conjecture on the lower bound for the permanent 
of doubly stochastic matrices in 1926 [2 which was later proved (in 1981) by 
Egoritchev and by Falikman as exposed by Knuth in [8]. 

Theorem 1 (Van der Waerden's Conjecture) Let A be a doubly stochastic n- 
square matrix, then 

per(A) ^ — 

n n 

with equality iff a% j — 1/n for all 1 ^ n. 
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Mine's Conjecture. There is also a result for the upper bound of the per- 
manent, due to an conjecture originally due to Mine [14]. The conjecture was 
first proved by Bregman in 1973 and a simpler proof is due to Schrijver |20j . 

Theorem 2 (Mine's Conjecture) Let A be a n-square matrix with values in 
{0,1} and non-zero sums Ti, 

n 

pei(A) < n( r -) 1/n 

i=i 

There are only few matrices for which an explicit formula for the permanent 
is available, the derangements being one instance. In fact [12] , 

n z r 

per(zl n + J) = rt! — 

r=0 r ' 

which, for the derangements, gives with z — — 1, 

per( J - /„) = „!(1 - I + I - ... + (-!)«!) 

Pcrmanents can be used to evaluate the number of Su-Doku grids. The scarce 
results known about permancnts, however, yield only information on upper 
bounds of the number of grids rather than their exact number which was essen- 
tially computed by brute force in [BJ. 

2.3 An implementation of Knuth's "Dancing Links" algo- 
rithm 

In a famous paper [?] Donald Knuth proposed an algorithm and a very efficient 
implementation for the exact cover problem. While the paper expands on 
its application to pentominoes, tetrasticks and to the queens problems, the 
algorithm itself, which Knuth called Algorithm X "for lack of a better name" , 
has a much broader scope. Through proper formulation of Su-Doku grid and 
puzzle problems, it proved efficient at enumerating grids and solving problems 
of various sizes. 

Knuth's first insight is to point that the matrix representation of the exact 
cover or exact hitting set problems makes it a good candidate for backtracking. 
Algorithm X is a simple expression of a generic backtrack process. 

Knuth's backtracking algorithm for the exact hitting set problem. 

Algorithm X is a nondeterministic algorithm, defined on a given matrix A of Os 
and Is. Citing from Knuth's paper: 

If is empty, the problem is solved; terminate successfully. 
Otherwise choose a column deterministically . 
Choose a row such that nondeterministically . 
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Include in the partial solution. 
For each such that , 

- delete column from matrix ; 

- for each such that , delete row from matrix . 
Repeat this algorithm recursively on the reduced matrix . 

The nondeterministic choice of a row means that all such rows are succes- 
sively (or in parallel) selected for inclusion into the partial solution, the algo- 
rithm proceeding essentially in an independent way on these rows. The choice 
of the column c, on the other hand impacts the execution time and exploration 
path of the algorithm. Any systematic rule for choosing a column in the pro- 
cedure will find all solutions. Certain rules, however, work better than others. 

In the Su-Doku experiments we studied two such rules: the random rule, 
where the column is chosen at random in the reduced matrix, and the shortest 
rule, where the column having the smaller number of Is is selected. While for 
enumeration tasks these options make no real difference, as Algorithm X in this 
case behaves basically as a trial and error procedure, we found that for puzzles, 
the shortest rule always outperformed the other one. This is also the result 
of experiments ran on another well-known combinatorial puzzle, the Langford's 
problem. 

Knuth's "Dancing Links" implementation. In the original paper, Knuth 
also proposed a very efficient implementation of Algorithm X based on doubly- 
linked circular lists. Each element in the matrix A is represented as a structured 
object with pointers to the previous and next elements in the same row (left and 
right), to the previous and next elements in the same column (up and down) and 
an extra-pointer to a column header structure which keeps track of the column 
name, its size (the number of Is) and additional metric information which can 
be useful to monitor the performance of the algorithm. 

typedef struct col { 

struct col *1; 

struct col *r; 

struct cell *u; 

struct cell *d; 

int size; 

char name [8] ; 

ClientDataPtr clientData; 
} *colPtr; 

typedef struct cell { 
struct cell *1; 
struct cell *r; 
struct cell *u; 
struct cell *d; 
colPtr c; 
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} *cellPtr; 



The I and r fields of the column headers link remaining columns in the 
reduced matrix which need to be covered. Global variables point to the circular 
list of columns and to the partial solution: 

static struct col S_Header; 
static cellPtr *S_Covering; 

With these data structures, the concrete implementation of Algorithm X is 
as follows: 

search( k ) : 

If S_Header.r == S_Header, print the current solution and return. 

Otherwise choose a column structure . 

Cover column . 

For each row in while , 

- set S_Covering[k]=; 

- for each in while , cover column ; 

- search ( k+1 ) ; 

- set =S_Covering[k] , and ; 

- for each in while , uncover column . 
Uncover column and return. 

The search procedure is initially called with k = to enumerate all solutions. 

Knuth's second insight is used to implement the cover/uncover function 
which are used to remove and reinstall columns in the matrix. Knuth observed 
that the "atomic" remove operations in a doubly-linked circular list: 

x.r.l — x.l andxJ.r = x.r 

are simply reversed, provided the x data structure is kept intact, by the subse- 
quent operations: 

x.r.l — x and x.l.r = x 

which will put back x in the circular list. 

The cover operation uses the first set of operations to remove a column first 
from the header list and then to remove all rows in c's own circular list from 
the other column lists they are in: 

int cover ( colPtr col ){ 
int updates = 0; 

/* Remove col from header list */ 
col->r->l = col->l; col->l->r = col->r; 
updates++; 

/* Remove all rows in col list from other col lists they are in */ 
cellPtr cell, rowCell; 
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for( cell = col->d; cell != (cellPtr) col ; cell = cell->d ){ 

for( rowCell = cell->r; rowCell != cell; rowCell = rowCell->r ){ 
rowCell->d->u = rowCell->u; rowCell->u->d = rowCell->d; 
rowCell->c->size -= 1; 
updates++ ; 

} 

} 

return updates; 

} 

The function also keeps track of counters for statistical purposes and decre- 
ments the column size in the header. The uncovering operation is symmetric, 
taking place in precisely the reverse order of the covering operation: 

/* 

* uncover - Inverse cover 
*/ 

void uncover ( colPtr col ){ 
/* Inserts all row cells */ 
cellPtr cell, rowCell; 

for( cell = col->u; cell != (cellPtr)col; cell = cell->u ){ 

for( rowCell = cell->l; rowCell != cell; rowCell= rowCell->l ){ 
rowCell->c->size += 1; 

rowCell->d->u = rowCell; rowCell->u->d = rowCell; 

} 

} 

/* Inserts in header list */ 
col->r->l = col; col->l->r = col; 

> 

The disconnected then reconnected links perform what Knuth called a "dance" 
which gave its name to this implementation known as the "Dancing Links" . 

The running time of the algorithm is essentially proportional to the number 
of times it applies the remove operation, counted here with the updates variable. 
It is possible to get good estimates of the running time on average by running 
the above procedure a few times and applying techniques described elsewhere 
by Knuth [?] and Hammersley and Morton [?] (so called "Poor Man's Monte 
Carlo"). 

Deriving cover matrices from representation matrices. We now turn 
to the proper formulation of Su-Doku questions for calculations by the Dancing 
Links algorithm. Considering an elementary constraint in the size 2 Su-Doku 
puzzle, for instance on one row of the 4x4 grid, its matrix representation as in 
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12. H is simply: 



M = 



V i 



1 \ 



1 / 



In M each column stands for a variable in the alldifferent constraint; there are 
four of them, one for each cell in the grid's row under consideration. The 
rows represent each of the possible four values 1, 2 ... 4 of all these variables. 
Remember that in a permutation constraint all variable domains are the same 
and their size is equal to the number of variables, here four. 

And of course, should we be interested in a single permutation constraint, the 
number of solutions, which as mentioned above is expressed by the permanent 
of this "all-Is" matrix, is simply, in this case, the number of permutations of 
four elements, i.e. 4! = 24. 

Now in order to obtain the A matrix for the Dancing Links algorithm, we 
augment the matrix M with the fact that each variable must have only one 
value. This is captured by four additional columns, one for each variable, 
containing a 1 for a given (row) value assigned to the variable and otherwise: 



X\ x 2 x 3 X4 C\ C 2 C 3 C 4 



1 1 

1 1 

1 1 

1 1 



= 3 



In this table the Xi column is the i-th variable in the C permutation constraint, 
the Cj column represents the i-th position in the constraint and each row a 
value from the shared domain {1, 2, 3, 4}. The A matrix has 8 columns and 16 
rows. 

More generally speaking, for a size n permutation constraint the A matrix 
counts n 2 rows and n + n — 2n columns. 

For a complete Su-Doku grid, there is one such constraint per line, per hie 
and per block. In addition, variables are shared between constaints, each one 
appearing in 3 constraints. Let us consider a Su-Doku of size n, which contains 
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n 4 , cells in n 2 lines by n 2 files, and n 2 blocks. The full size A matrix for the 
Dancing Links algorithm has n 4 + n 4 + ra 4 + n 4 = 4n 4 columns, one for each of 
the cells, and n 2 for each of the line, file and block in the grid. It also has n 6 
rows, one row for each of the n 2 possible value for each of the n 4 cells. The 
following table indicate the matrices sizes for different Su-Doku problems: 

n Matrix(rows x cols) 

2 64 x 64 

3 729 x 324 

4 4096 x 1024 

These sizes are small enough for the algorithm to perform satisfactorily on 
modern PCs. 

Enumerating grids and solving puzzles with the "Dancing Links". 

Having augmented the matrix to prepare it for the Dancing Links algorithms 
we are now ready to put the algorithm through different chores. 

In order to enumerate all Su-Doku valid grids we simply run the search ( 
) procedure with the appropriate A matrix as above. Of course, while the 288 
solutions of the size 2 Su-Doku grids are quickly enumerated, the size 3 grid takes 
evidently too long to list. Interestingly enough, size 2 variations of Su-Doku 
grids, such as diagonal Su-Doku grids where in addition one requires that all 
numbers in both diagonals to be different - adding two additional permutation 
constraints to the existing set, captured by 2n 2 additional columns in the A 
matrix - can also be enumerated by the same procedure. 

In order to solve puzzles, we need to remove from matrix A the rows cor- 
responding to the givens in the puzzle. In our implementation, this is simply 
another parameter file to the command line. If there arc k such givens in 
the puzzle, k rows are initially added to the partial solution and the procedure 
search ( k ) is called. The algorithm then proceeds, as above, to enumerate 
all solutions to the puzzle. It can be used to validate a puzzle, making sure 
that it has only one solution. 

2.4 Experimentation and results 

Enumerating size-2 Su-Doku grids. Running the Dancing Links algorithm 
on the 64 by 64 size-2 Su-Doku A matrix, produces the first of the 288 solutions 
almost immediately: 

Read 64 columns from sud2.mat 

Read 64 rows from file sud2.mat 

[16] New covering 1/1 in sees, usees: 



Depth Covers Backtracks Degrees 

37 1 4 

1 25 1 2 

2 22 1 2 

3 16 1 1 
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4 


28 


1 


3 


5 


16 


1 


1 


6 


19 


1 


2 


7 


10 


1 


1 


8 


10 


1 


1 


9 


16 


1 


2 


10 


10 


1 


1 


11 


7 


1 


1 


12 


16 


1 


2 


13 


10 


1 


1 


14 


7 


1 


1 


15 


7 


1 


1 


Total 


256 


16 





Estimation of solution path: 
7620 

The sud2.mat file is the A matrix for the size-2 Su-Doku grid. The trace 
table shows the depth, i.e. the value of k which indicates the depth in the 
backtrack tree; the cover count, which is the number of elementary remove 
operations in the circular lists; the number of backtracking steps at each depth 
level; and the degree, the number of children nodes explored at each level. 
Finally the estimation of the average number of operations to reach a solution 
is printed according to the "Poor Man's Monte Carlo" method. 

Counting Su-Doku grids. The algorithm can be used to count the number 
of Su-Doku grids, here for the sizc-2 grid: 



Read 


64 


columns from sud2.mat 




Read 


64 


rows 


from file sud2.mat 




1 




16 


7620 


7620 


2 




16 


7620 


15240 


3 




16 


5316 


20556 


4 




16 


5316 


25872 


5 




16 


7620 
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16 
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41112 


7 




16 
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48732 
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16 


7620 


56352 
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16 


5316 


61668 


10 




16 


5316 
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11 




16 
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74604 


12 




16 
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82224 
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16 
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16 
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15 




16 
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102780 


16 
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16 
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115716 


18 




16 
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24 


16 
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164448 


25 


16 


7620 


172068 
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16 
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179688 
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16 


5316 
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185004 


28 


16 


5316 


190320 


29 


16 


7620 


197940 


30 


16 


7620 
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205560 


31 


16 
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213180 


32 


16 
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220800 


33 
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5316 


226116 


34 


16 
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231432 


35 


16 
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16 
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16 
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46 


16 
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47 


16 
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323580 


48 


16 
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328896 


49 


16 
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334212 
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16 
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341832 


51 


16 
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349452 


52 


16 
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357072 


53 


16 
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54 


16 
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370008 


55 


16 
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56 


16 
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382944 


57 


16 
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390564 


58 


16 
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398184 


59 


16 
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405804 
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16 


5316 


411120 


61 


16 


5316 


416436 


62 


16 
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424056 


63 


16 
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431676 


64 


16 
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439296 
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16 
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16 
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70 


16 
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71 


16 
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72 


16 
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16 
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513900 


76 


16 
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77 


16 
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526836 


78 


16 
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534456 


79 


16 


7620 


542076 


80 


16 


7620 


549696 


81 


16 


5316 


555012 


82 


16 
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560328 


83 


16 


7620 
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567948 


84 


16 
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575568 
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86 
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89 
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16 
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95 


16 
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96 


16 
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97 


16 
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98 


16 
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99 


16 


5316 
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16 
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104 


16 


7620 


714144 


105 


16 


5316 


719460 


106 


16 


5316 


724776 


107 


16 


7620 


732396 


108 


16 


7620 


740016 


109 


16 


5316 


745332 


110 


16 


7620 


752952 



30 



Ill 


16 


7620 


760572 


112 


16 


7620 


768192 


113 


16 


7620 


775812 


114 


16 


5316 


•-70-1 -t fi 

781128 


115 


16 


5316 


"~70 f A A A 

786444 


116 


16 


7620 


794064 


117 


16 


7620 


801684 


118 


16 


7620 


809304 


119 


16 


7620 


816924 


120 


16 


5316 


822240 
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16 


5316 


1053672 


155 


16 


7620 


1061292 
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16 
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183 


16 


7620 


1253916 
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6852 

The first column prints the solution's number, the second the depth reached 
(always 16 for size-2 grid), the third one the estimation of the tree size, and the 
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last one the cumulated estimation used to provide the average on the last line. 
Here, in average, 6852 nodes are explored in the backtrack tree problem space. 



Further results. The Dancing Links implementation makes it very easy to 
experiment with several variations of the Su-Doku grids and puzzles. Adding 
further constraints to the problem is simply a matter of adding columns to the 
A matrix used by the algorithm. In the diagonal variant, for instance, where a 
Su-Doku grid is considered valid if, in addition, all numbers in both diagonals 
are also different, 2n columns are added to the A matrix to account for the n 
possible positions of each 1 . . . n 2 figure in each diagonal. 

In this variant, running the Dancing Links for enumeration yields the 48 
unique solutions for a size 2 diagonal Su-Doku problem (a 4-by-4 grid) and an 
average of 3666 nodes explored in the backtrack tree problem space. Note that 
the exploration space/time complexity is roughly halved on this instance. 

On a different track, the Dancing Links algorithm was successfully used for 
experimenting with the Langford problems which combinatorial nature, ulti- 
mately relying on permutation constraints, lends it perfectly to Dancing Links- 
based study. 
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